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Abstract 

The  effect  of  diverse  sources  of  uncertainties  and  the  intrinsically  multi-scale  nature  of  physical 
systems  poses  a  considerable  challenge  in  their  analysis.  Such  phenomena  are  particularly  critical  in 
material  systems  where  the  microstructural  variability  and  randomness  at  different  scales  have  a 
significant  impact  on  the  macroscopic  behavior  of  the  system.  Toward  this  goal,  during  the  period  of 
this  grant,  we  have  developed  a  sophisticated  though  efficient  and  accurate  multiscale  stochastic 
framework  for  uncertainty  quantification.  A  methodology  is  first  developed  to  incorporate 
topological  uncertainties  in  microstructures  using  a  non-linear  data-driven  model  reduction  technique. 
This  framework  seamlessly  allows  for  accessing  the  effects  of  microstructural  variability  on  the 
reliability  of  macro-scale  systems  and  provides  an  accurate  stochastic  input  model  into  our  stochastic 
system.  Next,  to  solve  the  resulted  stochastic  partial  different  equations  (SPDEs),  an  adaptive  sparse 
grid  collocation  technique  has  been  developed.  In  this  framework,  we  construct  the  stochastic 
collocation  points  based  on  the  function  being  represented,  thus  avoiding  computational  overhead. 
We  further  extended  this  framework  to  include  the  High  Dimensional  Model  Representation  (HDMR) 
technique  in  the  stochastic  space  to  represent  the  model  output  as  a  finite  hierarchical  correlated 
function  expansion  in  terms  of  the  stochastic  inputs  starting  from  lower-order  to  higher-order 
component  functions.  In  this  way,  we  can  address  the  stochastic  high  dimensional  problem  for  the 
first  time  in  this  area. 

We  applied  this  framework  for  the  design  of  general  materials  processes  under  uncertainty  including 
the  robust  design  of  deformation  processes  of  polycrystalline  materials.  We  developed  a  unique 
data-driven  strategy  to  encode  the  limited  information  on  initial  texture  and  grain  distribution  in 
deformation  processes  and  represent  it  in  a  finite-dimensional  framework.  We  have  developed  the 
methodology  to  produce  the  probabilistic  distribution  of  the  macro-scale  properties  of  the  material 
subjected  to  a  specific  process  induced  by  the  uncertainty  in  initial  microstructure  (texture  and  grain 
size  distribution). 


1  Summary  of  accomplishments 

Some  of  key  achievements  during  the  period  of  this  grant  are  given  below.  More  recent 
developments  (work  in  press  or  under  review  -  will  only  briefly  introduced). 

•  Development  of  a  non-linear  model  reduction  strategy  to  constmct  stochastic  input  models  of 
meso-scale  topology  variations  based  on  limited  data  (emphasis  on  polycrystalline  materials). 
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•  Development  of  an  adaptive  hierarchical  sparse  grid  collocation  algorithm  for  stochastic  partial 
differential  equations. 

•  Development  of  a  stochastic  multiscale  paradigm  to  address  simultaneously  the  effects  of 
randomness  and  multiscale  nature  of  physical  systems. 

•  Development  of  a  stochastic  optimization  technique  for  robust  design  of  deformation  processes  of 
polycrystalline  metals. 

•  Development  of  a  surrogate  stochastic  model  for  accelerating  multiscale  estimation  in  Bayesian 
inference  approaches. 

•  Development  of  a  maximum  entropy  approach  for  predicting  macroscopic  property  variability 
induced  by  uncertainty  in  initial  microstructure  in  deformation  processes. 

•  Development  of  an  HDMR  framework  for  representing  the  input/output  relation  of  complex 
systems  in  high-dimensions. 

Some  of  these  contributions  are  summarized  below  and  for  additional  details  the  relevant  references 
should  be  contacted. 

2  Data-driven  methodology  to  construct  stochastic  input  models  of  meso-scale 
topological/property  variations  [1][2][3] 

The  importance  of  performing  stochastic  analysis  on  heterogeneous  media  necessitates  the 
development  of  realistic  input  models  of  the  microstructural  features.  The  thermal,  mechanical  and 
chemical  behavior  of  microstructures  is  highly  anisotropic  and  heterogeneous,  depending  on  the 
randomness  of  features  of  importance.  For  instance,  orientation  of  the  crystals  as  well  as  the  nature  of 
the  grain  boundaries  represents  sources  of  randomness  in  polycrystalline  materials.  Knowledge  of  the 
topology/property  variation  of  say,  a  polycrystalline  material  is  usually  known  only  in  a  statistical 
sense  (in  terms  of  say,  grain  size  distribution  and  the  texture  map).  To  provide  reliable  failure  criteria 
for  critical  applications  involving  such  materials,  it  becomes  imperative  to  access  this  variability  in 
properties,  quantify  it  and  predict  its  effect  on  the  performance  of  the  system.  Stochastic  analysis  of 
random  heterogeneous  media  provides  information  of  significance  only  if  realistic  input  models  of 
the  topology  and  material  property  variations  are  used.  In  this  work,  we  have  introduced  a  framework 
to  construct  such  input  stochastic  models  for  the  topology,  thermal  diffusivity  and  permeability 
variations  in  heterogeneous  media  using  a  data-driven  strategy.  Given  a  set  of  microstructure 
realizations  (input  samples)  generated  from  given  statistical  information  about  the  medium  topology, 
the  framework  constructs  a  reduced-order  stochastic  representation  of  the  topology  and  material 
properties.  This  problem  of  constructing  a  low-dimensional  stochastic  representation  of  property 
variations  is  analogous  to  the  problem  of  manifold  learning  and  parametric  fitting  of  hyper-surfaces 
encountered  in  image  processing  and  psychology. 

2.1  Non-linear  model  reduction  strategies 

The  principal  component  analysis  (PCA)  based  model  reduction  scheme  constructs  the  closest  linear 
subspace  of  the  high-dimensional  input  space  [1].  This  is  a  fairly  good  approximation  when  the  number 
of  data  sets  is  limited.  But  as  the  amount  of  data  available  increases,  PCA  based  techniques  tend  to 
consistently  over-estimate  the  actual  dimensionality  of  the  space  [1].  This  is  primarily  due  to  the  fact 
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that  the  space  of  all  plausible  microstructure  distributions  (Q)  is  a  non-linear  space.  The  large  number  of 

available  data-points  effectively  populates  this  non-linear  space.  In  this  context,  non-linear 
transformation  strategies  offer  the  possibility  of  constructing  optimal  low-dimensional  representations 
of  this  space.  Here,  assume  that  M  plausible  microstructures  are  available,  each  of  which  is 

represented  as  a  high  dimensional  vector.  Using  the  unordered  data|kf|^  ,  kf  eQ,  the  problem  of 

interest  is  to  find  a  low-dimensional  parameterization  offi  ,  i.e.  a  set  r  e  K",  Af  <s  ,  such  that  there  is 
a  one-to-one  correspondence  between  Q  and  T . 

The  solution  strategy  is  based  on  the  so-called  principle  of  “manifold  learning”  [2].  The  basic 
strategy  is  to  show  that  this  set  of  unordered  points  lie  on  a  manifold  embedded  in  a  high-dimensional 
space.  That  is,  Q  is  a  manifold  embedded  in  a  high-dimensional  space.  The  mathematical  framework  is 
then  to  “unravel  and  smoothen”  this  manifold  and  represent  it  as  a  smooth  low-dimensional  curve  T . 
This  “unraveling  and  smoothing”  corresponds  to  a  topological  transformation  that  preserves  some 
notion  of  the  geometry  of  the  manifold.  The  framework  essentially  boils  down  to  two  mathematical 
steps; 

•  Defining  the  appropriate  manifold  on  which  this  high-dimensional  data  lie  on  and  identifying 
some  properties  of  the  manifold,  and 

•  Defining  the  appropriate  transformation  that  results  in  the  low-dimensional  equivalent  space. 

By  defining  an  appropriate  distance  function/) between  two  points  inD,  we  construct  a  metric  space 
{0.,D)  .  The  first  constraint  that  we  impose  during  the  construction  of  a  transformation  g :  Q  ^  r  is  that 

Q  is  topologically  well-behaved,  i.e.  it  is  smooth  and  has  no  holes.  This  can  be  ensured  by  showing 
thatO.  is  compact  [2].  The  next  step  is  to  choose  a  geometric  feature  of  the  manifold  and  construct  a 
transformation  that  keeps  this  feature  invariant  under  the  transformation.  The  key  notion  is  that  by 
keeping  specific  geometrical  features  of  the  embedded  manifold  invariant  one  can  construct  a  low¬ 
dimensional  representation  that  is  equivalent  to  the  manifold.  A  natural  choice  of  a  geometric  feature  is 
the  distance  metric.  This  results  in  an  isometric  mapping  to  transform/)  into T .  The  important  idea  is 
that  the  distance  that  encodes  the  geometric  information  about  the  non-linear  manifold  in  the  geodesic 
distance.  The  geodesic  distance  reflects  the  true  geometry  of  the  manifold  embedded  in  the  high¬ 
dimensional  space. 

Construction  of  T  reduces  to  finding  a  low-dimensional  representation  {Y,. }  of  the  given  data  points 
kf,...,k^  such  that {yJ is  isometric  tokf,...,k®  based  on  the  geodesic  distances  between  the  points. 
Denote  the  intrinsic  geodesic  distance  between  points  inQ  by  is  defined  by 

(kf  =  inf  {length  (/)}  (1) 

where  y"  varies  over  the  set  of  smooth  arcs  cormecting  kf  and  kf  .  The  essential  step  in  preserving 
distances  is  to  first  compute  the  pair-wise  distance  in  the  manifold  between  all  the  data  points 
|kf,...,k"  I .  This  brings  an  apparent  paradox  that  has  to  be  resolved:  The  intrinsic  geometry  of  the 

manifold  is  unknown.  To  detect  the  geometry  as  a  means  of  constructing  a  low-order  representation, 
we  utilize  the  notion  of  the  geodesic  distance.  But  the  construction  of  the  geodesic  distance  requires 
some  idea  of  the  underlying  geometry  of  the  manifold  (see  Eq.(l)).  An  approximation  of  the  geodesic 
distance  is  required  to  proceed  further.  Such  an  approximation  is  provided  via  the  concept  of  graph 
distance.  The  unknown  geodesic  distances  in  Q  between  the  data  points  are  computed  in  terms  of  a 
graph  distance  with  respect  to  a  neighborhood  graph  G  constructed  on  the  data  points.  This 
neighborhood  graph  G  is  very  straightforward  to  construct.  Two  points  share  an  edge  on  the  graph  if 
they  are  neighbors  (with  the  edge  length  being  proportional  to  the  distance,/)  between  them)  [2].  For 
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points  close  to  each  other,  the  geodesic  distance  is  well  approximated  by  the  distance  measured  by  the 
distance  metric  D  .  This  is  because  the  curve  can  be  locally  approximated  to  be  a  linear  patch,  and  the 
geodesic  distance  between  two  points  on  this  patch  is  the  straight  line  distance  between  them.  On  the 
other  hand,  for  points  positioned  faraway  from  each  other,  the  geodesic  distance  is  approximated  by 
adding  up  a  sequence  of  short  hops  between  neighboring  points.  These  hops  can  be  computed  easily 
from  the  neighborhood  graph  G  .  This  approximation  asymptotically  matches  the  actual  geodesic 
distance  (Eq.(l))  as  the  number  of  samples,  M  increases  [2]. 

The  computation  of  the  approximate  geodesic  distance  between  all  pairs  of  points  is  a  key  step  in 
the  framework.  By  appropriately  defining  the  distance  metric  D ,  we  have  encoded  all  the  information 
about  the  geometry  of  the  manifold  (utilizing  the  given  set  of  data  |kf  ,...,k"  | )  into  the  geodesic 

distance  matrix  (denoted  asM).  The  estimation  of  the  low-dimensional  representation  ofjkf  ,...,k®  | 
can  now  be  posed  as: 

Find  a  configuration  o//iomts'{Y,,...,Y^},  Y,.  such  that  these  points  yield  a  Euclidean  distance 
matrix  whose  elements  are  identical  to  the  elements  of  the  geodesic  distance  matrixM  .  That  is,  find 
{Y,  }“  such  that||Y,. The  principle  of  Multi-dimensional  scaling  (MDS)  can  subsequently  be 

used  to  compute  the  set  of  low-dimensional  points  that  best  represent  the  high-dimensional  points.  The 
MDS  procedure  essentially  computes  the  eigen-decomposition  of  the  geodesic  matrix  and  sets  the  low¬ 
dimensional  points  as  linear  combinations  of  the  largest  eigenvectors  of  the  geodesic  matrix. 

The  fact  that  Q.  is  compact  ensures  that  r  is  a  convex,  connected  region  inR"  [2].  This  provides  a 
natural,  elegant  way  of  constructing  T  from  the  M  low-dimensional  points  {Y,  }“  : 

r  =  {Y€R'"|YG  convex  hull{Y,,...,Y^}}  (2) 

The  intrinsic  dimensionality  N  of  the  low-dimensional  representation  can  be  estimated  by  using  a 
variant  of  the  Breadwood-Halton-Hammersley  [2]  theorem  (a  powerful  result  in  geometric  probability) 
where  it  is  linked  to  the  rate  of  convergence  of  the  length-functional  of  the  minimal  spanning  tree  of  the 
geodesic  distance  matrix  of  the  unordered  data  points  in  the  high-dimensional  space. 

The  overall  steps  of  the  procedure  are  summarized  in  Fig  1.  Note  that  here  the  surrogate  space  of 
convex  hull  is  mapped  to  a  unit  d  dimensional  hypercube  to  allow  interfacing  this  procedure  with 
sparse  grid  collocation  techniques  discussed  below.  In  such  collocation  methods,  the  sampling  points 
are  defined  on  a  hypercube.  These  collocation  methods  have  been  shown  to  be  efficient  in  interfacing 
with  deterministic  solvers  of  e.g.  deformation,  diffusion,  flow,  thermal,  etc.  in  random  media,  thus 
allowing  modeling  the  effect  of  microstructural  uncertainty  on  material  properties. 
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Reconstruction  for  arbitrary 
points  on  the  hypercube 


Unit  hypercube 


Convex  hull  of  the  points: 
The  low-dimensional  model 


Figure  1:  The  various  steps  in  data-driven  model  reduction  of  polycrystal  microstructures.  The  high¬ 
dimensional  microstructures  are  mapped  to  a  low-dimensional  region.  This  convex  region  is  mapped 
to  a  unit  hypercube.  Each  sampling  point  on  this  hypercube  corresponds  to  a  microstructure  that 
needs  to  be  reconstructed  using  the  given  data.  This  model  is  then  served  as  the  input  into  SPDEs 
and  solved  using  sparse  grid  collocation  discussed  below. 


3  Hierarchical  adaptive  framework  for  the  solution  of  stochastic  partial  differential  equations 
(SPDEs)  [4]  [5] 

The  above  generated  N-dimensional  representation  of  the  stochastic  fine-scale  material  property  is 
utilized  as  an  input  stochastic  model  for  the  solution  of  SPDEs.  We  utilize  an  adaptive  sparse  grid 
collocation  strategy  for  constructing  the  stochastic  solution.  We  briefly  describe  the  development  of  the 
adaptive  sparse  grid  collocation  strategy  here.  More  details  are  given  in  our  recent  work  in  [4] [5].  This 
technique  is  general  as  it  can  be  applied  for  creating  a  high-dimensional  interpolant  (in  the  stochastic 
space)  to  the  solution  of  any  stochastic  PDE  physical  system  using  only  the  deterministic  solver  as  a 
black  box  simulator. 

The  collocation  method  requires  only  repetitive  calls  to  an  existing  deterministic  solver  similar  to  the 
Monte  Carlo  method  and  performed  better  than  pre-existing  spectral  techniques  when  the  number  of 
random  dimensions  is  high.  However,  both  the  SSFEM  (stochastic  spectral  finite  elements)  and 
sparse  grid  collocation  methods  utilize  global  polynomials  in  the  stochastic  space.  In  the  presence  of 
steep  gradients  or  finite  discontinuities  in  the  stochastic  space,  these  methods  converge  very  slowly 
or  even  fail  to  converge.  Later  [5],  we  extended  this  method  to  an  adaptive  sparse  grid  collocation 
strategy  (ASGC)  using  piecewise  multi-linear  hierarchical  basis  functions  wherein  the  concept  of 
hierarchical  surplus  is  used  as  an  error  indicator  to  detect  regions  of  discontinuities  in  the  stochastic 
space.  The  basic  idea  here  is  to  use  a  piecewise  linear  hat  function  as  a  hierarchical  basis  function  by 
dilation  and  translation  on  equidistant  interpolation  nodes.  Then  the  stochastic  function  can  be 
represented  by  a  linear  combination  of  these  basis  functions.  The  corresponding  coefficients  are  just 
the  hierarchical  increments  between  two  successive  interpolation  levels  (hierarchical  surpluses).  The 
magnitude  of  the  hierarchical  surplus  reflects  the  local  regularity  of  the  function.  For  a  smooth 
function,  this  value  decreases  to  zero  quickly  with  increasing  interpolation  level.  On  the  other  hand, 
for  a  non-smooth  function,  a  singularity  is  indicated  by  the  magnitude  of  the  hierarchical  surplus.  The 
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larger  this  magnitude  is,  the  stronger  the  singularity.  Thus,  the  hierarchical  surplus  serves  as  a  natural 
error  indicator  for  the  sparse  grid  interpolation.  When  this  value  is  larger  than  a  predefined  threshold, 
we  simply  add  the  2N  neighboring  points  to  the  current  point.  A  key  motivation  towards  using  this 
framework  is  its  linear  scaling  with  dimensionality,  in  contrast  to  the  A-dimensional  tree  (2^)  scaling 
of  the  /?-type  adaptive  framework.  In  addition,  such  a  framework  guarantees  that  a  user-defined  error 
threshold  is  met.  In  particular  we  also  showed  that  it  was  rather  easy  with  this  approach  to  extract 
realizations,  higher-order  statistics,  and  the  probability  density  function  (PDF)  of  the  solution. 


3.1  Adaptive  sparse  grid  collocation  method 

The  basic  idea  of  this  method  is  to  utilize  a  finite  element  approximation  for  the  spatial  domain 
and  approximate  the  multi-dimensional  stochastic  space  using  interpolating  functions  defined  on  a 
set  of  collocation  points  eF  .  The  interpolation  can  be  constructed  by  using  either  full-tensor 

product  of  ID  interpolation  rule  or  the  so  called  sparse  grid  interpolation  method  based  on  the 
Smolyak  algorithm  [4].  Since  the  number  of  support  points  grows  very  quickly  as  the  number  of 
stochastic  dimensions  increases  in  the  full-tensor  product  case,  we  mainly  focus  on  the  sparse  grid 
method  and  discuss  the  proposed  adaptive  algorithm. 

3.1.1  Smolyak  algorithm 

The  Smolyak  algorithm  provides  a  way  to  construct  interpolation  functions  based  on  a  minimal 
number  of  points  in  multi-dimensional  space.  Using  Smolyak  method,  univariate  interpolation 
formulae  are  extended  to  the  multivariate  case  by  using  tensor  products  in  a  special  way.  This  provides 
an  interpolation  strategy  with  potentially  orders  of  magnitude  reduction  in  the  number  of  support  nodes 
required.  The  algorithm  provides  a  linear  combination  of  tensor  products  chosen  in  such  a  way  that  the 
interpolation  property  is  conserved  for  higher  dimensions. 

Let  us  consider  a  smooth  function  /  :[0,1]^^M  .  In  the  ID  case,  we  consider  the  following 

interpolation  formula  to  approximate  /  ;  =  with  the  set  of  support  nodes 

x‘={y;\y;  e  [0, 1]  for 7  =  1, 2, . . .  m,.  I ,  where  i  e  N,  a‘j  =  (f)  )  e  C ([0, 1])  are  the  interpolation  nodal  basis 

functions,  and  m,  is  the  number  of  elements  of  the  set  X‘ .  We  assume  that  a  sequence  of  the  ID  formula 
is  given  with  different  i .  In  the  multivariate  case  N  >1,  the  tensor  product  formulae  are 

j\  -1  Jn 

this  serves  as  building  blocks  for  the  Smolyak  algorithm. 

With  =  0,A'  =  U‘  -  U‘~\\ i  1=  (  +•••!„ for  i  e ,  and  q>N,  the  Smolyak  algorithm  constructs  the 
sparse  interpolant  „  as  [5] 

|i|<g  |i|=g 

' - V - ' 

To  compute  the  interpolant  4^  ^  from  scratch,  one  needs  to  compute  the  function  at  the  nodes  covered 
by  the  sparse  grid  ^  : 
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(5) 


The  construction  of  the  algorithm  allows  one  to  utilize  all  the  previous  results  generated  to  improve  the 
interpolation.  By  choosing  appropriate  points  for  interpolating  the  ID  function,  one  can  ensure  that  the 
sets  of  points  X'  are  nested  (X'  c  ) .  To  extend  the  interpolation  from  level  / -1  to  / ,  one  only  has  to  evaluate 

the  function  at  the  grid  points  that  are  unique  to  X'  ,  that  is,  atX^  =  X'  \X'“^  .  Thus,  to  go  from  an  order  q-\ 
interpolation  to  an  order  q  interpolation  in  N  dimensions,  one  only  needs  to  evaluate  the  function  at  the  differential 
nodes  XH^  ^  given  by 

|i|=9 

3.1.2  Choice  of  collocation  points  and  the  nodal  basis  functions 

The  advantage  of  choosing  equidistant  nodes  is  its  capability  for  allowing  adaptivity.  We  consider 
the  ID  interpolation  with  the  number  of  nodes  defined  asm,  =  1,  if  i  =  1;2'"‘  +1,  if  /  >  1 .  Then  the  supports 
nodes  are 


T'  = 


7-1 

m,  -1 


,  for 7  =  1, . . . ,  m, ,  if  m,  >  1;  0.5,  forj  =  1,  if  m,.  =  1 


(7) 


It  is  noted  that  the  resulting  grid  points  are  nested  and  it  is  called  “Newton-Cotes”  grid.  The  simplest 
choice  of  ID  basis  function  is  the  standard  linear  hat  function  [5].  a{Y)  =  l-\Y\,  if7e[-l,l];  0,  otherwise. 

This  mother  of  all  piecewise  linear  basis  functions  can  be  used  to  generate  arbitrary  a'  with  local 
support [f,'  +2‘“']by  dilation  and  translation,  i.e. 


a\  =  1  for  i  =  1,  and  a'.  = 


fl-(m,, -1)-\Y-Y;\,  if  \Y-Y;\<V{m,  -l) 
I  0,  otherwise 


(8) 


fori > land 7  =  l,...,m,. The  N-dimensional  multi-linear  nodal  basis  functions  can  be  constructed  using 
tensor  products  as  follows; 


(9) 


where  the  multi-index  j  =  e  N" and  j^,k  =  l,...,N ,  denotes  the  location  of  a  given  support  node 

in  the  k  -th  dimension. 


3.1.3  From  nodal  basis  to  multivariate  hierarchical  basis 

Let  us  consider  the  incremental  interpolant  formula  Eq.  (4).  This  formula  takes  advantage  of  the 
subset  property  of  the  nested  grid  points  X‘  a  X**' .  Here,  we  follow  closely  [5]  to  provide  a  clear 
development  of  the  derivation  of  the  hierarchical  basis  and  the  hierarchical  surpluses.  We  start  from  the 
ID  interpolating  formula  using  nodal  basis  as  discussed  in  the  previous  section.  We  have 
A'(/)  =  C'(/)-f/'“‘(/)  with  U‘{f)  =  Y,y,^^,a)-f[Y’),  andc/'-'(/)  =  C/'([/'-‘(/)),  we  obtain 


(10) 


Y‘gX‘ 


yUx^ 


Y‘gX‘ 


1 


(11) 


and,  since  f{Y‘j^-U‘  ‘ (/)(7^')  =  0, V7^'  e  X‘  ' ,  we  obtain 

A'(/)= 

recalling  that  X'^  =  X‘  \X‘“'  .  Clearly,  X'^  has  m‘^  =  w,  -m,_,  points,  since  X‘  a  .  For  simplifying  the 
notation,  we  consecutively  number  the  elements  in  and  denote  the  j  -th  point  ofx^  as  .  Then  we 
can  rewrite  the  above  equation  as 

m\ 

(12) 

7=1  ^ - V - ' 

Wj 

Here,  we  define  w/ as  the  ID  hierarchical  surpluses,  which  is  just  the  difference  between  the  function 
value  at  the  current  and  previous  interpolation  levels.  We  also  define  the  set  of  functions  a'  as  the 
hierarchical  basis  functions. 

For  the  multi-dimensional  case,  through  a  new  multi-index  set 

5,  eX'*  foij,  ,k  =  1,...,n},  (13) 

we  can  define  the  hierarchical  basis  as  jal^ :  j  e  <  i|. 

Now  we  apply  the  ID  Eq.  (12)  to  obtain  the  sparse  grid  interpolation  formula  for  the  multivariate 
case  in  a  hierarchical  form.  From  Eq.  (4),  we  obtain 

4-1, ).(/)=  Z 


For  smooth  functions,  the  hierarchical  surpluses  tend  to  zero  as  the  interpolation  level  increases.  On  the 
other  hand,  for  non-smooth  functions,  steep  gradients/fmite  discontinuities  are  indicated  by  the 
magnitude  of  the  hierarchical  surplus.  The  bigger  the  magnitude  is,  the  stronger  the  underlying 
discontinuity  is.  Therefore,  the  hierarchical  surplus  is  a  natural  candidate  for  error  control  and 
implementation  of  adaptivity. 

As  a  matter  of  notation,  the  interpolation  function  used  will  be  denoted  as  4^+^  ,  where  k  is  called 
the  level  of  the  Smolyak  interpolation.  We  consider  the  interpolation  error  in  the  space 

Fj^  •=(/  40,1]^  — continues,  OT;  <2,Vij,  (15) 

where  rnsNo  andD^^'is  the  usual  N -variate  partial  derivative  of  order  |in| :  •••07™") . 

Then  the  order  of  the  interpolation  error  in  the  maximum  norm  is  given  by  [4]  [5] 

l|/-4).(/)|L=o(^"|log2Mr-‘>),  (16) 

whereM  =  dini(/f(^,A())  is  the  number  of  interpolation  points. 

3.1.4  From  hierarchical  interpolation  to  hierarchical  integration 

Any  function  M  e  F  can  now  be  approximated  by  the  following  reduced  form  from  Eq  (14)  ; 
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(17) 


«(x,Y)  =  XX>v:(x)-a:(Y) 

|i|^?  je5i 

It  is  just  a  simple  weighted  sum  of  the  value  of  the  basis  functions  for  all  collocation  points  in  the 
current  sparse  grids.  Therefore,  we  can  easily  extract  the  useful  statistics  of  the  solution  from  it.  The 
mean  of  the  random  solution  can  be  evaluated  as  follows; 

E[M(x)]  =  XZ^i«-J«|(YyY  (18) 

|i|<g  jeSj  r 

where  the  probability  density  function  p{Y)  is  1  since  we  have  assumed  uniform  variable  variables  on  a 
unit  hypercube  [0,1]"  .  The  ID  version  of  the  integral  in  the  equation  above  can  be  evaluated 
analytically; 

a'j  (Y)dY  =  1,  if ;  =  1;  1  /  4,  if  i  =  2;  2'“' ,  otherwise.  (19) 

This  is  independent  of  the  location  of  the  interpolation  point  and  only  depends  on  the  interpolation 
level  in  each  stochastic  dimension  due  to  the  translation  and  dilation  of  the  basis  function.  Since  the 
random  variables  are  assumed  independent  of  each  other,  the  value  of  the  multi-dimensional  integral  is 
simply  the  product  of  the  ID  integrals.  Denoting  j^ai(Y>/Y  =  /i  ,  we  can  write  Eq.  (18)  as 

E[m(x)]  =  •  Thus,  the  mean  is  just  an  arithmetic  sum  of  the  hierarchical  surpluses  and  the 

|i|<^  je5j 

integral  weights  at  each  interpolation  points. 

To  obtain  the  variance  of  the  solution,  we  need  to  first  obtain  an  approximate  expression  forw^ ,  i.e. 
(x,  Y)  =  ^  vi  (x)  •  a\ (Y) .  Then  the  variance  of  the  solution  can  be  computed  as; 

|i|^?  je5i 


Var [m  (x)]  =  e[m"  (x)]  - (e[m  (x)])'  =  X  S  W '  “ 


(20) 


3.1.5  Adaptive  sparse  grid  interpolation 

In  this  section,  we  will  develop  an  adaptive  sparse  grid  stochastic  collocation  algorithm  based  on  the 
error  control  of  the  hierarchical  surpluses.  Before  discussing  the  algorithm,  let  us  first  introduce  some 
notation.  The  ID  equidistant  points  of  the  sparse  grid  in  Eq.  (7)  can  be  considered  as  a  tree-like  data 
structure  as  shown  in  Figure  2.  It  is  noted  that  special  treatment  is  needed  here  from  level  2  to  level  3. 
For  the  boundary  nodes  in  level  2,  we  only  add  one  point  along  the  dimension  (there  is  only  one  son 
here  instead  of  two  sons  for  all  other  levels  of  interpolation).  Then  we  can  consider  the  interpolation 
level  of  a  grid  point  Fas  the  depth  of  the  tree  Z)(F) .  For  example,  the  level  of  a  point  0.25  is  3.  Denote 
the  father  of  a  grid  point  asE’(F) ,  where  the  father  of  the  root  0.5  is  itself,  i.e.  E’(0.5)  =  0.5.  We  denote 
the  sons  of  a  grid  point  Y  =  (}],...  F„)  by 

Sons(Y)  =  {s  =  (5„5„...,5^)|(E(5,),5„...,5^)  =  Y,...,  or  (5„5,,...,E(5^))  =  y}  (21) 

From  this  definition,  it  is  noted  that,  in  general,  for  each  grid  point  there  are  two  sons  in  each 
dimension,  therefore,  for  a  grid  point  in  a  A  -dimensional  stochastic  space,  there  are  2N  sons.  It  is  also 
noted  that,  the  sons  are  also  the  neighbor  points  of  the  father.  Recall  from  the  definition  of  grid  points 
from  Eq.  (7)  and  the  definition  of  hierarchical  basis  that  the  neighbor  points  are  just  the  support  nodes 
of  the  hierarchical  basis  functions  in  the  next  interpolation  level.  By  adding  the  neighbor  points,  we 
actually  add  the  support  nodes  from  the  next  interpolation  level,  i.e.,  we  perform  interpolation  from 
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level  I  i  I  to  level  |  i  |  +1  .  Therefore,  in  this  way,  we  refine  the  grid  locally  while  not  violating  the 
developments  of  the  Smolyak  algorithm  Eq.  (14)  .  The  basic  idea  here  is  to  use  hierarchical  surpluses 
as  an  error  indicator  to  detect  the  smoothness  of  the  solution  and  refine  the  hierarchical  basis  functions 
ai  whose  magnitude  of  the  hierarchical  surplus  satisfies  |wi|  >  5  .  If  this  criterion  is  satisfied,  we  simply 

add  the  2N  neighbor  points  of  the  current  point  from  Eq.  (21)  to  the  sparse  grid. 


0.5 


0.25 

/  \ 

0.125  0.375 

/\  /\ 


\ 

0.75 

/  \ 

0.625  0.875 

/\  /  \ 


Figure  2:  ID-tree  like  structure  of  sparse  grid 

Therefore,  let£>0be  the  parameters  for  the  adaptive  refinement  threshold.  We  propose  to  use  the 
following  iterative  refinement  algorithm  beginning  with  a  coarsest  possible  sparse  gridA^^  ^,  i.e.,  with 

the  N  -dimensional  multi-index  i  =  (l, . . . l) ,  which  is  just  the  point  (0.5, . . . ,  0.5) . 

( 1 )  Set  level  of  Smolyak  construction  k  =  Q . 

(2)  Construct  the  first  level  adaptive  sparse  grid  „  . 

•  Calculate  the  function  value  at  the  point(0.5,...0.5) . 

•  Generate  the  IN  neighbor  points  and  add  them  to  the  active  index  set. 

•  Set  k  =  k  +  \ . 

(3)  While  k  <  k^^  and  the  active  index  set  is  not  empty; 

•  Copy  the  points  in  the  active  index  set  to  an  old  index  set  and  clear  the  active  index  set. 

•  Calculate  in  parallel  the  hierarchical  surplus  of  each  point  in  the  old  index  set  according  to 

Here,  we  use  all  of  the  existing  collocation  points  in  the  current  adaptive  sparse  grid  ^ .  This 
allows  us  to  evaluate  the  surplus  for  each  point  from  the  old  index  set  in  parallel. 

•  For  each  point  in  the  old  index  set,  if  |  wi  |  >  ^  ; 

Generate  2N  neighbor  points  of  the  current  active  point  according  to  Eq  (2 1 ) ; 

Add  them  to  the  active  index  set. 

•  Add  the  points  in  the  old  index  set  to  the  existing  adaptive  sparse  grid  „  .  Now  the 
adaptive  sparse  grid  becomes  A„^^  „ . 

•  k  =  k  +  l 

(4)  Calculate  the  mean  and  the  variance,  the  PDF  and  if  needed  realizations  of  the  solution. 


Figure  3  shows  some  representative  results  in  this  area. 
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Figure  3:  The  figure  depicts  the  standard  deviation  of  temperature  profile  across  a  random  Ag-W 
composite  microstructure  due  to  its  topological  uncertainty. 


4  An  adaptive  high  dimensional  stochastic  model  representation  technique  for  the  solution  of 
stochastic  partial  differential  equations  [6] 

Although  the  current  adaptive  collocation  problem  can  solve  a  wide  range  of  problems,  there  are  still 
some  difficulties  when  addressing  the  problem  with  random  heterogeneous  media.  As  is  well  know, 
in  realistic  random  heterogeneous  media  often  we  deal  with  a  very  small  correlation  length  and  this 
result  in  a  rather  high-dimensional  stochastic  space  with  nearly  the  same  weights  along  each 
dimension.  In  this  case,  all  the  current  methods  are  not  applicable.  Toward  this  end,  in  [6],  a 
computational  methodology  was  developed  to  address  the  solution  of  high-dimensional  stochastic 
problems.  It  utilizes  high-dimensional  model  representation  (HDMR)  technique  in  the  stochastic 
space  to  represent  the  model  output  as  a  finite  hierarchical  correlated  function  expansion  in  terms  of 
the  stochastic  inputs  starting  from  lower-order  to  higher-order  component  functions.  HDMR  is 
efficient  at  capturing  the  high-dimensional  input-output  relationship  such  that  the  behavior  for  many 
physical  systems  can  be  modeled  to  good  accuracy  only  by  the  first  few  lower-order  terms.  An 
adaptive  version  of  HDMR  was  also  developed  to  automatically  detect  the  important  dimensions  and 
construct  higher-order  terms  using  only  the  important  dimensions. 

HDMR  represents  a  function  in  high-dimensions  in  the  following  form: 

/=i  (22} 

l</i  <•  •  -<4  <N 
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Here  is  the  zeroth-order  component  function  which  is  a  constant  denoting  the  mean  effect.  The 
first-order  component  function  is  a  univariate  function  which  represents  individual 

contribution  to  the  output  /(T)  .  It  is  noted  that  (f^)  is  general  a  nonlinear  function.  The  second- 
order  component  function is  a  bivariate  function  which  describes  the  interactive  effects 

of  variables  Y.^  andl^^  acting  together  upon  the  output  /(Y)  .  The  higher-order  terms  reflect  the 

cooperative  effects  of  increasing  number  of  input  variables  acting  together  to  impact  f  The  last  term 
gives  any  residual  dependence  of  all  input  variables  cooperatively  locked  together  to  affect  the  output 
/(Y)  .  Once  all  the  component  functions  are  suitably  determined,  then  the  HDMR  can  be  used  as  a 
computationally  efficient  reduced-order  model  for  evaluating  the  output.  This  is  the  same  idea  as  the 
stochastic  collocation  method  where  we  also  obtain  an  approximate  representation  of  /  (  Y)  . 

In  this  work,  the  CUT-HDMR  is  adopted  to  construct  the  response  surface  of  the  stochastic  solution. 
With  this  method,  a  reference  point  Y  =  (^l^,T2,---,kAf)is  introduced.  The  component  functions  of 
CUT-HDMR  are  explicitly  given  as  follows 

/o=/U).  y:(>')=/(Y)U^-/o 

The  basic  conjecture  underlying  HDMR  is  that  the  component  functions  arising  in  typical  physical 
problems  will  not  likely  exhibit  high-order  cooperativity  among  the  input  variables  such  that  the 
significant  terms  in  the  HDMR  expansion  are  only  those  of  low  order.  Therefore,  it  is  expected  that 
the  HDMR  expansion  will  converge  very  fast.  For  most  well-defined  physical  systems,  the  first-  and 
second-order  expansion  terms  are  expected  to  have  most  of  the  impact  upon  the  output  and  the 
contribution  of  higher-order  terms  would  be  insignificant. 

Within  the  framework  of  CUT-HDMR,  we  can  write  it  in  a  more  general  form  as 
/(Y)=  I/.(Y.)=  2;Z(-ir"/(Y.U,,  (24) 

uc£)  ucD  vcu 

for  a  given  setu  c  D,  where  Z)  :=  |l, . . . ,  A^}  denotes  the  set  of  coordinate  indices  and  /0  ( Y^  )  =  f^. 
Here  Y„  denotes  the  |  u  |  -dimensional  vector  containing  those  components  of  Y  whose  indices 
belong  to  the  set  u ,  where  |  u  |  is  the  cardinality  of  the  corresponding  set  u . 

Therefore,  the  Y-dimensional  stochastic  problem  is  transformed  to  several  lower-order  |  v  |  - 
dimensional  problems  /  ( Y^  )y=y\y  which  can  be  easily  solved  by  the  ASGC  as  introduced  in  the  last 
section: 

/(Y)=ES{-ir"  Z  Z^-oKY.)  (25) 

ucD  vcu  ||i||<ZV-i-^  j 

where  1 1  i  ||=  I’l  H - 1-/^,  w'^are  the  hierarchical  surpluses  for  different  sub-problems  indexed  by  v  and 

al(Y^)  is  only  a  function  of  the  coordinates  which  belong  to  the  set  v.  It  is  noted  that  the 
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interpolation  level  q  may  be  different  for  eaeh  sub-problem  according  to  their  regularity  along  the 
particular  dimensions  which  is  controlled  by  the  error  threshold  £• . 


Interpolation  is  done  quickly  here  through  simple  weighted  sum  of  the  basis  functions  and  the 
corresponding  hierarchical  surpluses.  In  addition,  it  is  also  easy  to  extract  statistics  as  introduced  in 
ASGC  by  integrating  directly  the  interpolating  basis  functions.  Let  us  denote 


vcu 


(26) 


as  the  mean  of  the  component  function  .  Then  the  mean  of  the  HDMR  expansion  is  simply 
E[/(Y)]  =  .  To  obtain  the  variance  of  the  solution,  we  can  similarly  construct  an 

approximation  forn^  and  use  the  formula  Far[M(x)J  =  (x)J-^E[M(x)]j  . 


We  have  considered  two  issues  related  to  adaptivity.  At  first  the  component  functions  are  computed 
using  the  adaptive  sparse  grid  collocation  method  (ASGC)  [5].  The  error  in  ASGC  is  controlled  by 
the  user  based  on  the  values  of  the  hierarchical  surpluses  and  hierarchical  basis  functions.  By 
integrating  HDMR  and  ASGC,  it  is  computationally  possible  to  construct  a  low-dimensional 
stochastic  reduced-order  model  of  the  high-dimensional  stochastic  problem  and  easily  perform 
various  statistic  analysis  on  the  output.  The  second  level  of  adaptivity  is  to  decide  on  the  fly  which 
component  functions  to  compute.  Note  that  for  high-dimensional  problems  even  the  computation  of 
all  two-body  terms  is  computationally  very  expensive. 


At  first,  we  try  to  find  the  important  dimensions.  To  this  end,  we  always  construct  the  zeroth-  and 
first-order  HDMR  expansion  where  the  computational  cost  is  affordable  even  for  very  high- 
dimensions.  In  this  case,  a  weight  is  defined  as: 


7,  = 


IKi-iIm., 


(27) 


where  ^  f^{Y^)dY.  and  the  norm  is  defined  in  the  spatial  domain  when  the  output  is  a 

function  of  spatial  coordinates.  Then  we  define  the  important  dimensions  as  those  whose  weights  are 
larger  than  a  predefined  error  threshold  .  Now  the  set  D  only  contains  these  important  dimensions 


instead  of  all  the  dimensions.  However,  not  all  the  possible  terms  are  computed.  Instead,  we 
adaptively  construct  higher-order  component  functions  increasingly  from  lower-order  to  high-order 
in  order  to  reduce  the  computational  cost  in  the  following  way.  For  each  computed  higher-order  term 
/„ ,  I  u  |>  2 ,  a  weight  is  also  defined  as 


ikJ 


Aw 


vg5,|v|<|u|-1 

L^{D) 

(28) 


It  measures  the  relative  important  with  respect  to  the  sum  of  current  integral  value  which  has  already 
been  computed  in  set  S  from  previous  order.  Similarly,  the  important  component  functions  are 
defined  as  those  whose  weights  are  larger  than  the  predefined  error  threshold  6^ .  We  put  all  the 
important  dimensions  and  higher-order  terms  in  to  a  set  T ,  which  is  called  the  important  set.  When 
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adaptively  constructing  HDMR  for  each  new  order,  we  only  calculate  the  term  whose  indices 
satisfy  the  admissible  relation: 

ueD and  v c u  => v e T  (29) 

In  other  words,  among  all  the  possible  indices,  we  only  want  to  find  the  terms  which  can  be 
computed  using  the  previous  known  important  component  functions.  In  this  way,  we  find  those  terms 
which  may  have  significant  contribution  to  the  overall  expansion  while  ignoring  other  trivial  terms 
thus  reducing  the  computational  cost  for  high-dimensional  problems. 

In  this  work,  we  examined  physical  processes  (hydrodynamic  transport,  deformation,  etc.)  in  random 
heterogeneous  media  and  have  reported  examples  of  up  to  500  random  dimensions  (note  this  is  the 
highest  stochastic  dimension  problem  that  is  currently  reported  in  the  literature  based  on  non  Monte 
Carlo  based  approaches).  Figures  4,5  provide  some  typical  results  for  flow  in  random  heterogeneous 
media. 


Figure  4:  The  left  figure  shows  the  standard  deviation  of  the  v-velocity  component  across  y =0.5. The 
right  figure  shows  the  convergence  of  PDF  at  one  point.  The  problem  here  refers  to  flow  (squared 
domain)  in  random  media  using  and  exponential  kernel  for  the  log-permeability  (with  high 
variability  defined  by <7^  -2.0 and  500  stochastic  dimensions.  The  parameter  O^is  used  to  control 
the  adaptive  selection  of  the  critical  dimensions. 


4  Decoupled  stochastic  multiscale  framework  [7] 

Based  on  all  the  previous  developed  methodologies  for  stochastic  problem,  we  are  able  to  apply  them 
to  various  applications  related  to  random  materials.  First,  we  developed  a  stochastic  variational 
multiscale  formulation  to  incorporate  uncertainties  in  multiscale  material  systems.  In  this  scheme,  a 
stochastic  analogue  to  the  mixed  multiscale  finite  element  framework  is  used  to  formulate  the 
physical  stochastic  multiscale  process.  For  the  effective  resolution  of  the  multiscale  problem,  the 
solution  was  split  using  an  additive  decomposition  into  its  coarse  scale  and  fine  scale  parts.  We 
employ  the  local  conservation  assumption  through  which  we  convert  the  global  sub-grid  problem 
into  a  set  of  local  sub-grid  problems.  This  is  the  first  time  that  a  multiscale  variational  technique  has 
been  applied  for  stochastic  PDFs. 

In  [7],  we  applied  this  framework  to  analyze  flow  through  random  heterogeneous  media  when  only 
limited  statistics  about  the  permeability  variation  are  given.  Linear  and  non-linear  model  reduction 
techniques  are  used  to  convert  the  limited  information  available  about  the  permeability  variation  into 
a  viable  stochastic  input  model.  An  adaptive  sparse  grid  collocation  strategy  is  used  to  efficiently 
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solve  the  resulting  stochastic  partial  differential  equations.  However,  our  mathematical  developments 
in  this  context  are  very  generic  in  nature  and  can  be  easily  extended  to  other  applications. 

4.1  Problem  definition 

Denote  the  domain  as  Da  R"’" ,  where  is  the  number  of  space  dimensions.  The  characteristic 

length  scale  ofD  isL .  Denote  the  length  scale  of  permeability  fluctuation  as  / .  In  the  problems  that  we 
are  interested  in  solving,  the  characteristic  length  of  the  domain  is  a  couple  of  orders  of  magnitude 
larger  than  the  characteristic  length  scale  of  the  permeability  fluctuations /«:  i .  We  are  interested  in 
evaluating  the  pressure,  pand  velocity,  u  in  the  domain,  D .  The  variables (u(x),/?(x)) depend  on  the 

(multiscale)  permeability  distribution,  k{x)  in  the  domain.  However,  the  complete  permeability 

distribution  is  unknown.  Only  some  limited  statistics  and/or  snapshots  of  the  permeability  are  given. 
This  limited  information  available  to  characterize  the  permeability  necessitates  assuming  that  the 
permeability  is  a  realization  of  a  random  field.  This  is  mathematically  stated  as  follows: 

Letfibe  the  space  of  all  allowable  permeability  variations.  This  is  our  event  space.  Every  point 
k  =  {k{TL  ,®),  Vx  G  A®  e  o}  in  this  space  is  equiprobable.  Consequently,  we  can  define  a  a  -algebra C  and 

a  corresponding  probability  measureP:C^[0,l]to  construct  a  complete  probability  space (Q,C,yc>) of 

allowable  permeability.  To  make  this  abstract  description  amenable  to  numerical  simulation,  a  finite 
dimensional  approximation/representation  of  this  abstract  set  is  necessary.  Various  data-driven 
strategies  to  represent  the  setD  are  discussed  in  Section  4.  The  stochastic  permeability  is  represented  as 
k{'s.,co)  =  k{'s.,Y^,...Yf^)  =  k{x,\) ,  where T , •  •  •  are  uncorrelated  random  variables. 

The  pressure  and  velocity  are  characterized  by  the  following  set  of  equations 

V-u(x,Y)  =  /(x) 

^  ^  '  (30) 

u(x,Y)  =  -/:(¥)  V/i(x,Y) 

Here,  the  source/sink  term / (x)  is  taken  to  be  deterministic. 

The  basic  idea  is  to  solve  the  problem  on  a  coarse  spatial  discretization  while  taking  into 
account  the  fine-scale  variation  in  the  stochastic  permeability.  In  the  next  section,  we  detail  a  stochastic 
extension  to  the  variation  multiscale  method  to  solve  this  problem.  The  stochastic  multiscale 
formulation  is  based  on  the  multiscale  formulation  detailed  in  the  paper  [7]. 


4.2  Variational  multiscale  formulation 

For  the  problem  to  be  physically  relevant,  we  assume  that  the  stochastic  permeability  A:  is  positive 
and  uniformly  coercive.  As  stated  in  Section  2,  the  abstract  representation  ofA: (•,<»)  in  Q  is  replaced  by  a 

more  tractable  finite  dimensional  representation  A: (•,¥)  ,  withFGrcR"  .  Corresponding  to  the 
probability  measure  P:C^[0,l]  ,  we  denote  the  equivalent  probability  measure  yo:r ->^[0,1]  .  The 
governing  equations  for  the  velocity  and  pressure  given  in  the  mixed  form  are  as  follows: 


k  'u-l- V/7  =  0 
V-u  =  / 


(31) 


with  the  following  boundary  conditions  p  =  p„  ondD^  andu  n  =  ondD^ .  Without  loss  of  generality,  we 

assume  that  the  boundary  conditions  are  deterministic  and  that  the  Neumann  condition  is  homogeneous 
=  0  on  dD^ . 

The  next  step  is  to  introduce  the  appropriate  function  spaces  in  which  the  velocity  and  pressure  lie. 
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We  introduce  the  following  tensor  function  spaces  W  =  5” Off  =  with  the  inner  product 

defined  as{p,p)  =  \pf^  :=  p[\)^ ^p^ dx =  S ® H  =  ll (T)® H [div,D) .  We  will  also  use  the  function 
space  defined  as  V  =  {u ;  u  e H,u(-,Y)  n  =  0  on  VY  e  r} .  Thus,  the  problem  can  be  written  in  mixed 
variational  form;  Find(u,jti)  c  Vx  W  such  that 


(v,A:''u)-(V-v,/7)  =  -(v-n,/7^),  Vv  e  V 
(w,V-u)  =  (w,/),  VweW 


(32) 


where(/,g)  is  defined  as  J^(/yC'(Y)|^^  fgdx  . 

In  the  variational  multiscale  approach,  the  exact  solution  u  is  assumed  to  be  made  up  of  contributions 
from  two  different  (spatial)  scales  namely,  the  coarse-scale  solutionu^(x,  )that  can  be  resolved  using  a 

coarse  (spatial)  mesh  and  a  sub-grid  solution (x,  )  such  that:  u  =  +u^ and p  =  p^+  Pf  This  additive 

sum  decomposition  induces  a  similar  decomposition  for  the  spatial  part  of  the  fine-scale  tensor-product 
function  spaces  into  a  direct  sum  of  a  coarse-scale  and  a  sub-grid  tensor-product  function  spaces,  e.g. 
W  =  W^©  .  The  main  idea  is  to  develop  models  for  characterizing  the  effect  of  the  sub-grid  solution 

u^(x,  )on  the  coarse  scale  solution  and  to  subsequently  derive  a  modified  coarse  scale  formulation  that 
only  involves  (x,  •) .  The  additive  decomposition  provides  a  way  of  splitting  the  fine-scale  problem 

given  by  Eq.  (32)  into  a  coarse-scale  problem  and  a  sub-scale  problem.  Testing  against  the  coarse-scale 
test  functions  results  in  the  coarse-scale  variational  problem;  Find(u^,/7j  e  V^x  such  that 

( V. ,  («.  +  «/ ))  -  ( V  •  V, ,  ))  =  - ( •  n,  J ,  Vv,  € 

(33) 

(w^,V-(u^+u^))  =  (w^,/),  Vw^gW^ 

Similarly  testing  against  the  sub-scale  test  functions  results  in  the  sub-scale  variational  problem:  Find 
(u^,/>^)g  V^x  such  that 


( v^ ,  A:-'  (u,  +  ))  -  ( V  •  v^ ,  ))  =  - ( v^  •  n,  ) ,  Vv^  G 

(34) 

(w,,V-(u,+u^))  =  (w^,/),  Vw^gW, 

The  key  is  to  solve  Eq.  (34)  for  and  construct  a  functional  representation  of  the  sub-scale  variation, 
u^and p^m  terms  of  the  coarse-scale  variation, =0(uj,  ='F(uJ  .  This  representation  can  be 
subsequently  used  to  remove  explicit  dependence  of  and  p^  in  Eq.  (33)  as: 

(v^,r'(u^+(l)(uJ))-(V-v^,(/7^+'P(uJ))  =  -(v^-n,;?„),  Vv^gV^ 
(w„V-(u,+(l)(uJ))=(w,,/),  Vw,gW, 

The  key  problem  is  now  to  solve  Eq.  (34)  over  each  coarse  scale  element  and  utilize  this  sub-grid 
stochastic  solution  to  solve  the  stochastic  coarse  scale  equation.  For  a  detailed  discussion  on  the 
solution  of  the  two  scale  stochastic  problems  please  refer  to  the  paper  [7]. 


4.3  The  stochastic  multiscale  framework 

The  abstract  framework  to  solve  the  stochastic  multiscale  problem  defined  by  Eq.  (30)  is  as 
follows: 
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variational  slop^iastt^  ^  /  Adaptive  sparse  grtd^ 
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framework 
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Stochastic  multiscale  velocity  and 
pressure 


Figure  5:  Schematic  of  the  developed  stochastic  multiscale  framework 

Fig.  6  shows  the  effect  of  uncertainty  in  multiscale  permeability  on  the  flow  through  a  random 
porous  medium. 


Figure  6:  Left:  The  multiscale  log-permeability  distribution  in  the  domain.  Right:  Mean  contour  of 
the  stochastic  coarse  scale  x-direction  flux. 


5  Predicting  property  variability  of  polycrystals  induced  by  microstructural  uncertainty:  A 
maximum  entropy  approach  [8] 

The  quantification  and  propagation  of  uncertainty  in  process  conditions  and  initial  microstructure  on 
the  final  product  properties  in  a  deformation  process  were  investigated.  The  stochastic  deformation 
problem  was  modeled  using  the  sparse  grid  collocation  approach.  The  ability  of  the  method  in 
estimating  the  statistics  of  the  macro-scale  microstructure-sensitive  properties  and  constructing  the 
convex  hull  of  these  properties  is  shown  through  examples  featuring  randomness  in  initial  texture 
and  process  parameters.  A  data-driven  model  reduction  methodology  together  with  a  maximum 
entropy  approach  is  used  for  representing  randomness  in  initial  texture  in  Rodrigues  space. 
Comparisons  are  made  with  the  results  obtained  from  the  Monte-Carlo  method.  In  modeling  the 
texture  evolution,  the  random  initial  texture  was  represented  as  a  random  field.  The  available 
information  on  initial  microstructure  provided  as  a  set  of  x-ray  diffraction  images  is  rarely  enough  to 
completely  define  the  aforementioned  random  field.  In  this  situation,  one  needs  to  resort  to  the 
maximum  entropy  approach  in  which  the  random  field  is  constructed  such  that  the  entropy  of  the 
information  it  conveys  is  maximized.  The  method  used  in  this  work  is  fairly  general  and  as  the 
known  information  on  the  microstructure  increases  it  can  be  easily  incorporated  in  approximating  the 
random  field. 
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5.1  Constitutive  problem  and  texture  evolution 

Consider  a  point  in  the  reference  fundamental  region  that  corresponds  to  a  particular  crystal 
orientation.  In  an  appropriate  kinematic  framework,  the  total  deformation  gradient  is  decomposed 
into  plastic  and  elastic  parts,  F  =  FT’’ ,  where  F®  is  the  elastic  deformation  gradient  and  F’’ ,  the 
plastic  deformation  gradient,  with  detF^  =  1 .  The  constitutive  relation  is  given  by 

f  =  Z"[F"]  (36) 

where  T  is  the  second  Piola-Kirchhoff  stress  tensor,  F  is  the  fourth-order  anisotropic  elasticity 
tensor  expressed  in  terms  of  the  crystal  stiffness  parameters  and  the  orientation  r  and 

=  ^{f^  F”"  -  /j  .  The  re-orientation  velocity  is  found  as  follows: 

V  =  ^  =  ^[a  +  {cd-r^r  +  dj^r^  (37) 

where  ris  the  orientation  (Rodrigues’  parameterization)  and  0  represents  the  spin  vector  defined  as 
a  =  vect^F^'F*'^) ,  where  F^'is  evaluated  through  the  polar  decomposition  of  the  elastic  deformation 

gradient  F"  as  F"  =F"C/F 

Consider  a  macroscopic  material  point  and  an  associated  underlying  microstructure  M  discretized  by 
a  finite  element  grid.  Each  point  on  this  underlying  grid  corresponds  to  a  different  crystal  orientation 
F  .  At  each  point  on  the  grid,  the  crystal  lattice  frame  e^  is  related  to  the  sample  reference  frame  e^  by 

e^  =  Re^ .  Due  to  crystal  symmetry,  the  orientation  F  is  not  unique.  Restricting  the  Rodriguez  domain 

to  a  fundamental  zone  that  reflects  the  crystal  symmetry  leads  to  a  one  to  one  correspondence 
between  the  points  on  the  Rodriguez  space  and  the  crystal  orientation. 

The  Rodrigues-Frank  axis-angle  parameterization  is  used  as  a  convenient  scheme  to  represent  F  .  The 
parameterization  is  derived  from  the  natural  invariants  ofF  :  the  axis  of  rotations  and  the  angle  of 
rotation^ .  The  angel-axis  parameterization,  r ,  is  obtained  by  scaling  the  axis  n  by  a  function  of  the 

angles  asr  =  «  /($■).  In  the  particular  case  of  Rodrigues’  paramterization,  the  function  is  defined  as 
/(f)  =  tanf^  . 

The  Lagrangian  scheme  for  the  ODF  evolution  is  used.  The  evolution  of  the  ODF  is  governed  by  the 
ODF  conservation  equation  and  is  given  in  the  Lagrangian  form  as  follows 

^ -I- V  •v(5',t)  =  0  (38) 

where  v(5',t)  is  the  Lagrangian  re-orientation  velocity  of  the  crystals  and  the  Lagrangian  form  of  the 
ODF,  ,  is  subjected  to  J^j',0)  =  Jq  (x)  as  the  initial  condition. 

See  Figures  6  for  the  schematic  of  the  problem. 
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Figure  7:  Schematic  view  of  the  effects  of  uncertainty  in  initial  texture  on  the  final  material 
properties.  The  error  bars  on  the  effective  stress/strain  response  at  a  material  point  shown  are  due  to 
the  uncertainty  in  initial  texture  and  variability  in  processing,  (bottom).  While  these  calculations  are 
at  a  material  point  of  a  polycrystal,  they  pave  the  way  for  computing  the  property  variability  in  a 
workpiece  during  processing  induced  from  lack  of  information  on  the  microstructure  of  the  initial 
workpiece. 


5.2  Problem  deflnitions:  Process  and  texture  uncertainty 

Consider  a  complete  probability  space  (Q,F,P)  where  Q  is  the  event  space,  F  the  cr -algebra,  and 
P  ;  F  ^  [0,  1]  is  the  probability  measure.  The  uncertainty  in  the  problem  we  consider  comes  from:  (a) 
variation  in  the  velocity  gradient  representing  the  variation  in  process  parameters:  e  Q  and 

(b)  variation  in  the  initial  texture:  ^  (.s,  <y) ,  x  g  91,  ft)  g  Q  . 


The  velocity  gradient  is  written  in  terms  of  various  deformation  modes  such  as  tension/compression, 
plain  strain  compression,  shear  and  rotation.  The  coefficients  of  these  terms  can  be 

assumed  as  random  variables  to  represent  variation  in  process  conditions. 
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The  incompressibility  condition  is  assumed  here  and  only  eight  component  of  L  are  independent  and 
hence  P  consists  as  well  of  eight  components. 
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One  can  use  a  random  field,  ^  ( j",  <y) :  ^  e  91,  <»  e  Q  to  represent  the  variability  of  the  initial  texture. 
The  stochastic  partial  differential  equation  for  the  evolution  of  texture,  A{^s,t,co) 
:91x[0,r]xQ^M^  -u{0}  ,  can  be  written  such  that  for  P-almost  everywhere 

dA(s,t,a))  . 

— ^ - -  +  A{^s,t,a>)V  ■v[s,t,a))  =  0  (39) 

In  this  work,  we  used  a  maximum  entropy  (MaxEnt)  principle  to  seek  a  joint  probability  distribution 
of  the  random  texture. 


5.4  Probability  distribution  of  the  random  variables  using  the  maximum  entropy  (MaxEnt) 
principle 

Let 7  =  (}^,...,Fjv)be  the  set  of  random  variables  for  which  the  probability  density  function is 

unknown.  This  probability  density  function  is  assumed  as  a  map  from  D  a  to  =  [0,  +oo[  where 
D  is  the  support  of  Py  and  is  defined  previously  as  the  convex  hull 
(  Z)  =  7  e  convexhull  cz  R}  )  of  all  admissible  values  of  7|,...,7^  .  Any  probability 

function  should  satisfy  the  following  constraints  in  order  to  be  acceptable  for  this  problem. 

E(f{Y))  =  M  (40) 

where  M  =  , . . . ,  )  is  a  given  vector  in  with  h  being  the  number  of  constraints  defined  in 

£(J'(®))  =  0 

E{Y,{m)Yj{o,))  =  S,j_  i,J  =  l,...,N 

E  is  the  expectation  and  7^/(7)  is  a  given  measurable  mapping  from  R^  to  R*  .  These 

conditions  are  the  result  of  specific  properties  of  the  Karhunen-Loeve  expansion.  Hence  they  can  be 
written  as 

/(y)=(r,e(r)) 

M=(0„,c) 

where /i  =  +  A^(A^  + 1) / 2  ,  0^  =  (0,...,0)  e  R^  ,  e(7)  is  a  vector  in  formed  from  the 

diagonal  and  upper  triangular  part  of  the  matrix  77^  and  eis  a  vector  in  that  has  0  and  I’s 

as  its  elements. 


(41) 


(42) 


Let  us  define  the  information  entropy  iS'(/7)  of  the  probability  density  function  p  as 

5(;i)  =  -|;i(7)log(/?(7))J7  (43) 

if  S  is  the  convex  collection  of  all  probability  density  functions  defined  on  D  which  satisfy  the 
above  constraints  and  -qo|  is  nonempty,  then  the  Maximum  Entropy  principle 

consists  of  finding  the  probability  distribution  p  that  maximizes  the  information  entropy: 

Py  =  arg  max  Sip)  (44) 

The  MaxEnt  problem  can  be  posed  as  an  unconstrained  optimization  problem  using  Lagrange 
multipliers.  In  this  method,  the  constraints  are  incorporated  into  the  cost  function  as 
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C  =  S(p)  +  {X,E(f{Y))-M)^  (45) 

where  X  =  represents  the  Lagrange  multipliers.  Maximizing  this  cost  function  is 

equivalent  to  maximizing  the  entropy  and  satisfying  the  constraints.  Since  5  (/>)  is  a  concave 

functional  the  set  |P:  >  — coj  is  convex  and  uniqueness  of  PySiS  the  solution  of  Eq.(44)  is 

guaranteed.  It  can  be  shown  that  the  solution  of  the  above  problem  can  be  represented  as 

p{Y)  =  Z  exp(-(/l,/(r)}  (7)  V7  e  D  (46) 

where  Z  is  a  normalization  constant  and  I^(F)the  indicator  function;  I^(F)  =  1  if  FeZ)  and 
I^(F)  =  0,  otherwise.  The  Lagrange  multipliers  are  chosen  such  that  they  satisfy  the  constraints. 

when  the  number  of  constraints  is  small,  the  Lagrange  parameters  can  be  obtained  by  a  simple 
gradient  method  but  as  the  number  of  constraints  becomes  significant  a  dual  approach  can  be  used  in 
which  the  problem  is  posed  as  an  optimization  problem  in  terms  of  the  Lagrangian  parameters.  The 
dual  optimization  problem  can  be  written  as 

X*  =  argmin^(L) 

n 

where  2^  =  |^  exp /  (2))^  j  dY .  The  function satisfies  the  following  properties 

^  =  i  =  (48) 

OAi 

where  M.  is  defined  in  Eq.(42).  From  these  equations,  it  is  clear  that  the  solution  of  Eq.(47)  satisfies 

the  constraints  posed  in  Eq.  (40).  The  solution  of  the  constrained  optimization  problem  posed  by  the 
Maximum  Entropy  approach  has  the  parametric  form  shown  by  Eq.  (46)  where  X  can  be  inferred  by 

minimizing  the  dual  function  ^(L)  .  This  means  that  the  solution  of  the  dual  problem 
corresponds  to  the  p(a)  that  maximizes  the  entropy. 

After  obtaining  the  joint  probability  distribution  for  the  input  model  using  the  MaxEnt  approach,  then 
we  can  solve  the  problem  using  ASGC  or  HDMR  method  presented  before.  Through  this  method,  the 
convex  hull  of  properties  was  obtained  from  a  material  subjected  to  uncertain  process  parameters  and 
initial  texture.  This  can  be  important  for  providing  us  with  the  means  to  quantify  how  well  process 
conditions  and  microstructure  need  to  be  known  to  attain  desired  properties  but  also  to  identify  risks 
(e.g.  failure  probabilities)  affiliated  with  critical  values  of  the  material  properties.  To  the  best  of  our 
knowledge  this  is  the  first  time  these  concepts  have  been  explored  in  the  analysis  and  design  of 
polycrystalline  materials.  Figure  7  shows  a  convex  hull  of  Bulk  modulus,  shear  modulus  and  Young 
modulus  for  an  FCC  polycrystal  material.  In  essence  we  have  computed  all  feasible  properties  that 
one  should  expect  in  a  given  deformation  process  (tension,  compression,  etc.)  when  the  initial 
microstructure  is  random. 
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Figure  8:  The  convex  hull  of  Bulk  modulus,  shear  modulus  and  Young  modulus  for  an  FCC 
polycrystal  obtained  in  tension  for  random  initial  texture  (uncertainty  driven  by  data).  Extremal 
properties  can  be  identified  together  with  the  affiliated  probabilities.  These  unique  ideas  are  very 
important  not  only  for  design  under  uncertainty  but  also  for  failure  prediction  from  extremal 
scenarios. 


6  Microstructure  model  reduction  and  uncertainty  quantiflcation  in  multiscale  deformation 
processes  [9] 

We  developed  what  we  think  is  the  first  framework  for  stochastic  multiscale  deformation  processes. 
Including  the  underlying  microstructure  and  its  evolution  for  every  integration  point  on  macroscale  is 
essential  in  quantifying  the  effect  of  deformation  process  on  macroscale  properties.  A  reduced-order 
model  for  representing  the  data-driven  stochastic  microstructure  input  is  developed.  The  multiscale 
random  field  representing  the  random  microstructure  is  decomposed  into  few  modes  in  different 
scales  (the  Rodrigues  space  for  representing  texture  on  mesoscale  and  the  continuum  macroscale 
space).  Realizations  from  a  stochastic  simulation  are  used  to  obtain  a  small  number  of  modes 
approximating  the  stochastic  filed.  Then  a  bi-orthogonal  expansion  is  used  to  describe  the  variability 
of  the  initial  microstructure.  The  coefficients  of  the  polynomial  chaos  terms  in  this  expansion  are 
obtained  using  projections  of  the  random  modes  on  the  chaos  polynomials. 

Each  integration  point  on  the  macro  scale  is  associated  it  with  a  random  microstructure.  We  consider 
simultaneously  model  reduction  on  both  scales.  To  reduce  the  stochastic  dimensionality,  we  use  all 
microstructure  data  at  all  points  in  the  continuum  in  a  bio-orthogonal  KLE  expansion  that  allows  us 
to  reduce  the  number  of  random  variables  that  drive  the  multiscale  simulation.  In  essence  the  model 
accounts  for  the  microstructure  correlation  from  point  to  point  in  the  continuum  and  at  the  same  time 
produces  spatial  eigenfunctions  (similar  to  the  POD  model  reduction  of  PDEs). 


6.1  A  multiscale  reduced-order  model  of  the  uncertain  initial  microstructure 

This  section  provides  a  framework  to  obtain  a  reduced-order  model  for  the  underlying  random 
microstructure  field  (here  Orientation  density  function  (ODE)  defining  texture).  Assume  an 
random  field  a  (x, s,  co)  defined  on  a  probability  space  (Q,  F,  P)  . 

a(x,5,(y):Z)x91xQ— >]R  (49) 
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where  D  is  the  spatial  domain,  tR  is  the  fundamental  part  of  Rodriguez  space,  Q  is  the  set  of 
elementary  events  and  e  Q  is  the  vector  of  random  inputs.  One  can  use  the  Karhunen-Loeve 

expansion  to  express  this  field  by  a  bi-orthogonal  representation  in  the  form 

00 

a{^x,s,co)  =  a(x,  j')  +  a(x,  =  a  (x,x)  +  X,  ft))  (50) 

i=\ 

wherea  is  defined  as  a  (x,x)  =  ^d(x,  j',ft)))  and  (•)  is  the  averaging  operation  defined  below,  yO,.are 
eigenvalues  of  the  eigenvalue  problem  defined  later  on,  the  y/.  are  modes  strongly  orthogonal  in 
Rodrigues  space,  O,.  are  spatial  modes  weakly  orthogonal  in  space  with  respect  to  an  inner  product 
defined  as 

{f,g)-=\{f,g)dx  (51) 

D 

{f,g)  =  \f{(o)g{o))p{o))d(o  (52) 

where  /)(fti)  is  the  probability  distribution.  The  strong  orthogonality  of  modes  in  Rodrigues  space 
can  be  written  as 

and  the  weak  orthogonality  of  spatial  modes  can  be  written  as 

=  (53) 


By  minimizing  the  distance  (based  on  the  norm  defined  in  Eq.  (51) )  between  the  Karhunen-Loeve 
expansion  and  the  random  field,  one  ends  up  with 

=  (54) 

and  from  the  orthogonality  condition 


O 


If- 

(x,ft))  =  — ^1  d{x,s,o))y/^{^s^ds 
ylPi 


(55) 

(56) 

(57) 


Eqs.  (53)  and  (55)  lead  to  the  following  eigenvalue  problem 

PiWi{s)  =  lc{s,s')i//.{s')ds' 

where  the  covariance  C  is  defined  as 

c(x,x')=— [\ )  di,„  I  di,  I 

j=\  ,_=i  i„=i 

where  |  J.  \  is  the  Jacobian  determinant  of  the  element  4  ,  is  the  integration  weight  associated 
with  the  integration  pointy  ,  i®  number  of  integration  points  in  each  element,  is  the  number 
of  realizations,  n^i  is  the  number  of  elements  in  macroscale  and  a  is  a  column  vector  with  elements 
corresponding  to  integration  points  in  Rodrigues  space  and  x.^  represents  global  coordinate  of  the 
integration  pointy  in  macroscale. 


The  ODE  representing  the  texture  takes  positive  values.  Hence,  the  Karhunen-Loeve  expansion 
should  provide  us  with  positive  values.  To  obtain  a  positive  random  field,  one  can  use  the  Karhunen- 
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Loeve  expansion  for  the  fl(x, 5, ry)  =  log^^(x, assuming  ^4 (x, 5, ro)  >  >  0  almost 

surely.  The  proeess  A  can  be  reconstructed  as 


+  exp  (a  (x,  X,  ®))  =  +  exp  a  (x,  .s)  +  ^  ^t/^.  (s)  O.  (x,  co) 


(58) 


In  practice,  a(x,  j',<y)  :=  where  are  a  set  of  finite  number  of 

random  variables  and  refers  to  the  number  of  random  variables  considered  in  the  problem. 


(59) 


Next,  the  polynomial  chaos  decomposition  of  O.  (x,  <»)  can  be  written  as 

<D,  (x, (o)  :=  (D,  (x,g,{(o),..., (^y))  =  X 4  (^) (^) 

j 

where  the;/,  (ft))  =  7.  are  in  a  one-to-one  correspondence  with  the  Hermite  polynomials  in 

Gaussian  variables,  ^(ry)  is  the  vector  consisting  of  n,,  independent  Gaussian  random  variables 
(^^j , . . . ,  ^  and  the  coefficients  (p.j  (x)  can  be  obtained  from 

((l),(x,^)7.) 


4W=- 


(60) 


After  obtaining  the  reduced  model  on  the  two  scales  simultaneously,  ASGC  or  HDMR  can  be  used  to 
solve  the  resulted  SPDEs,  This  is  the  first  time  to  introduce  a  new  framework  that  makes  the 
otherwise  intractable  task  of  quantifying  the  effect  of  random  initial  texture  in  a  multiscale  problem 
feasible.  Figure  8  shows  a  representative  result  of  this  method. 
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Figure  9:  Left:  Schematic  view  of  the  multiscale  problem.  Each  point  on  this  underlying  coarse  grid 
corresponds  to  a  different  crystal  orientation  R  The  Rodrigues-Frank  axis-angle  parameterization  is 
used  as  a  convenient  scheme  to  represent  R  Model  reduction  is  conducted  simultaneously  on  the 
Rodrigues  space  and  the  coarse  grid.  Right:  Mean  and  variance  of  the  shear  modulus  obtained  from 
the  reduced-order  representation  of  texture.This  reduced  order  representation  is  then  used  in  a 
stochastic  multiscale  simulation  to  compute  the  variability  of  properties  in  the  final  product  (e.g. 
forged  product)  induced  from  the  uncertainty  in  the  initial  microstructure  [9J. 


24 


In  closing,  we  note  that  in  [13,14]  and  in  an  a  number  of  additional  forthcoming  publications,  we 
have  extended  all  of  the  stochastic  multiscale  polycrystal  material  models  discussed  above  to  include 
in  addition  to  texture  uncertainty  also  grain  size  uncertainty.  These  models  are  based  on  enhanced 
physical  modeling  to  account  for  the  effect  of  grain  size  distribution  on  macroscale  properties.  The 
model  reduction  techniques  (manifold  learning)  have  been  extended  to  work  on  high-dimensional 
microstructures  defined  in  both  the  texture  and  grain  size  random  space.  This  not  only  resulted  in 
appropriate  definition  of  distance  metrics  and  properties  of  the  microstructure  manifold  but  also  in  a 
number  of  necessary  extensions  of  the  MaxEnt  techniques  needed  to  compute  the  probabilistic 
distribution  of  the  underlying  random  variables.  At  the  end  of  this  work,  we  have  managed  to 
produce  the  convex  hull  of  all  material  properties  that  should  be  expected  (with  the  corresponding 
probabilities)  in  the  presence  of  microstructure  uncertainty.  This  has  implications  in  the  design  of 
new  materials,  on  prediction  of  rare  events  (e.g.  failure  due  to  extremal  properties),  etc. 
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New  Discoveries  and  Software 

(a)  Develop  an  HDMR  framework  to  address  the  high-dimensionality  of  stochastic  PDE  systems,  (b) 
Non-linear  reduced  order  model  that  could  capture  correlations  in  non-linear  spaces  and  efficiently 
represent/process  information  of  complex  stmctures,  (c)  developed  a  hierarchical  adaptive  sparse- 
grid  collocation  scheme  that  captures  the  cmcial  stochastic  dimensions  and  thus  solve  problems 
which  were  earlier  infeasible,  (d)  developed  a  variational  stochastic  multiscale  framework  for 
material  systems,  (e)  developed  a  non-intmsive  (collocation)  framework  for  design  of  complex 
systems  under  uncertainty  and  applied  it  to  the  design  of  deformation  processes  of  polycrystalline 
materials,  (f)  used  the  adaptive  sparse  grid  collocation  solver  as  a  surrogate  model  for  accelerating 
multiscale  Bayesian  inference  approaches  and  finally  (g)  developed  a  maximum  entropy  based 
framework  for  predicting  the  effects  of  uncertainty  in  initial  texture  on  macroscopic  property 
variability  in  deformation  processes  of  polycrystalline  materials. 

Several  of  the  computational  tools  developed  have  become  available  to  collaborators  at  other 
universities  as  well  as  to  government  laboratories. 
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